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ABSTRACT 


Different finite-element models previously set up for thermal analysis of the space shuttle orbiter structure were 
discussed and their shortcomings were identified. Element density criteria were established for the finite-element 
thermal modelings of space shuttle orbiter-type large, hypersonic aircraft structures. These criteria were based on 
rigorous studies on solution accuracies using different finite-element models having different element densities set up 
for one cell of the orbiter wing. Also, a method for optimization of the transient thermal analysis computer central 
processing unit (CPU) time was discussed. Based on the newly established element density criteria, the orbiter 
wing midspan segment was modeled for the examination of thermal analysis solution accuracies and the extent of 
computation CPU time requirements. The results showed that the distributions of the structural temperatures and 
the thermal stresses obtained from this wing segment model were satisfactory and the computation CPU time was 
at the acceptable level. The studies offered the hope that modeling the large, hypersonic aircraft structures using 
high-density elements for transient thermal analysis was possible if a CPU optimization technique was used. 

INTRODUCTION 

During the early stage of space shuttle development, the heat transfer analysis of the complex orbiter structure 
was carried out by using the so-called “plug” method. Only selected small, local regions of the orbiter were modeled 
with three-dimensional finite-difference (or lumped parameters) plug thermal models (fig. 1). Structural tempera- 
tures calculated from the plug models were then interpolated to obtain structural temperatures in the unmodeled 
regions. Approximately 90 percent of the temperature inputs used in the early days of orbiter thermal stress analysis 
was obtained by this type of interpolation. The procedure was tedious, laborious, and expensive. For structures with 
steep temperature gradients near the heat sinks (for example, the orbiter wing spar and nb caps), the structural tem- 
perature distribution obtained from the interpolation method could be erroneous and could cause inaccurate thermal 

stress predictions* 

After the development of finite-element heat transfer analysis computer codes such as the structural performance 
and resizing (SPAR) finite-element thermal analysis computer program (ref. 1), it became possible to model larger 
regions of the orbiter structure including the whole wing. The use of the finite-element method instead of the con- 
ventional lumped-parameter (or finite-difference) method enables the use of the same thermal model as a structural 
model for thermal stress calculations by simply removing elements set up for the thermal protection system (TPS) 
(not a major mechanical load-carrying structural component). In the past several years, Ko and others (refs. 2 to 
9) conducted extensive heat transfer and thermal stress analyses of the space shuttle orbiter using a series of finite- 
element models set up for three wing segments, one fuselage cross section, and the whole wing. These thermal 
models were used to calculate orbiter structural temperatures which were correlated with the actual flight-measured 
data during the initial orbiter tests of the space shuttle Columbia (refs. 2 to 8). The earlier thermal models set up for 
the orbiter structure were by no means perfect and had shortcomings in the light of element density, element size 
distributions, and extent of the region modeled. 

In the finite-element heat transfer and thermal stress analysis of small structural components, the element density 
may be freely increased at will to obtain highly accurate solutions without the worry that the limit of the computer 
memory core space might be reached or exceeded. However, in the finite-element thermal modeling of large, hy- 
personic aircraft structures (such as the space shuttle and space plane), the use of high-density finite elements could 
increase the number of radiation view factors tremendously, and could require prohibitive computer time and/or 
computer core space requirements in the transient heat transfer analysis using these radiation view factors. In the 
wing box-type structures, each time the number of radiation elements is doubled, the number of radiation view fac- 
tors would be nearly quadrupled. Thus, in the finite-element heat transfer analysis of large aerospace structures, 
the highest desirable element density is governed by the time requirement and/or by the memory capacity of the 
computer used. The corresponding structural model for thermal stress analysis using the same element density as 
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the thermal model has far less nodal points because the non-load-carrying heat shields are removed. It is, therefore, 
not the model to be used to set up the criteria for maximum allowable element density. Preliminary studies in this 
area were carried out by Ko and others (refs. 10 to 12) to investigate solution accuracies obtained from different 
finite-element models having different element densities set up for one wing ceil of the orbiter. These studies formed 
the foundation of the stale of the art in finite-element modelings of large, hypersonic aircraft structures. 

In this report, several past finite-element models set up for the orbiter structure are reviewed and their short- 
comings are identified. The criteria for the element density and element size distribution required for finite -element 
modeling of larger orbiter wing regions are discussed in detail. Finally, the report shows how to optimize the SPAR 
transient heat transfer analysis computation CPU time when the high-density elements are used. 

NOMENCLATURE 


c 

capacitance matrix 


C21 

two-node forced convection element 


C41 

four-node forced convection element 


CPU 

central processing unit 

- 

CQUAD2 

quadrilateral membrane and bending element 


CROD 

two-node tension-compression-torsion element 

— 

E23 

bar elements 


E25 

zero-length element for elastically connected, geometrically coincident joints 


E31 

triangular membrane element 

- 

E41 

quadrilateral membrane element 


E44 

quadrilateral shear panel element 

- 


radiation view factor from element i to element j 


FRSI 

flexible felt reusable surface insulation 

— 

H 

convection load vector 


HRSI 

high temperature reusable surface insulation 

- 

i 

integer, 1, 2, 3, 


JLOCS 

joint locations 

. = 

j 

integer, 1, 2, 3, ... 


K 

system matrix = Kk + K r + Kh 


K h 

convection matrix 


K k 

conduction matrix 

= 

K r 

radiation matrix 

— 

K21 

two-node line conduction clement 


K31 

three-node area conduction element 

«=■ 

K41 

four-node area conduction element 


K61 

six-node volume conduction element 

=- 

K81 

eight-node volume conduction clement 
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KTIME 

LRSI 

NASTRAN 

Q 

R 

R21 

R31 

R41 

RCONV 

RNITER 

RTIME 

SIP 

SPAR 

STS-5 

T 

T 2 

TPS 

Y 0 

['] 


control command in SPAR computer program to specify the time interval at which the system 
matrix if is to be factored 
low temperature reusable surface insulation 
NASA structural analysis 
source load vector 
radiation load vector 
two-node area radiation element 
three-node area radiation element 
four-node area radiation element 

control command in SPAR computer program to set the convergence criterion for the radiation load 
vector computations 

control command in SPAR computer program to set the maximum number of iterations permitted 
during the computations of radiation load vectors 
control command in SPAR computer program to specify the time increments at which the 
radiation load vectors are to be computed 

strain isolation pad 
structural performance and resizing 
space transportation system, flight 5 
absolute temperature 

final time of SPAR transient thermal analysis 

thermal protection system 

station in y-axis 

time derivative 

chordwise stress 

spanwise stress 

shear stress 


REVIEW OF PAST MODELS 


Wing and Fuselage Cross Sections 

Figures 2 and 3 show two typical past thermal models (WS240, FS877) set up for the orbiter wing midspan 
and midfuselage cross sections (ref. 5). The corresponding structural models were obtained by removing the TPS 
elements. In both figures, the sizes of the thermal and the corresponding structural models are compared. Notice 
that the structural models arc always simpler (less nodal points). These thermal models were extensively used to 
calculate orbiter structural temperatures which were correlated with the actual flight data obtained from initial orbital 
tests of the space shuttle Columbia (refs. 2 to 8). The WS240 model (fig. 2) has a reasonable number of elements 
to give satisfactory structural temperature distributions in the chordwise direction, but not in the spanwise direction. 
Also, the FS877 model (fig. 3), which is two-dimensional, can give good structural temperature distribution in the 
fuselage circumferential direction, but not in the fuselage axial direction. 
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Both of these thermal models, when converted to structural models by removing TPS elements, are not capable of 
yielding accurate thermal stress predictions because of the following reasons. It is well known that the magnitudes 
of the computed thermal stresses are sensitive to the modeling parameters (such as clement density, element size 
distributions, and extent of the region modeled) (refs. 1 1 and 12). For the case of biaxial stress fields which occur 
in orbiter-type structures (especially wing and fuselage skins), the manner in which the boundary conditions arc 
applied can greatly affect the magnitudes of calculated thermal stresses in the proximity of the traded boundaries or 
free edges. Because the thermal stresses attenuate to zero at free edges, the structural models must have dimensions 
sufficiently large enough to be beyond the influence of the edge effects. Both models mentioned earlier lack the 
capacity to fulfill those structural modeling requirements. 

Whole Wing 

In order to eliminate the modeling shortcomings that occurred in the WS240 models (thermal and structural), 
the whole wing of the orbitcr was modeled (refs. 8 and 9). Figure 4 shows the past thermal model, WING, set up for 
the orbitcr whole wing. In the figure the sizes of both the thermaTmodel and the structural model (TPS landing gear 
and wheel well door elements removed) are compared. In this case the number of nodes for the structural model is 
approximately 10 percent of that of the thermal model. 

Because of the fear of encountering and dealing with a tremendous numberof radiation view factors, the clement 
sizes used in the WING model were relatively coarse. The WING thermal model required 48,034 radiation view 
factor computations and can only give “roof’-shaped structural temperature distributions. Therefore, the thermal 
stresses calculated based on this type of temperature distribution may not give accurate predictions. Thus, finer 
elements are required for the whole wing model, for which the optimum manageable element density has to be 
determined first. 

ELEMENT DENSITY 

Because of the need of finding the optimum clement density and element size distribution for modeling orbiter 
whole wing-type structures, Ko and others (refs. 1 1 and 12) set up five different finite-element models (A, B, C, D, 
and E) with different element densities for one orbiter wing cell (located at midspan bay 3 (fig. 5)), and examined 
the effect of element density on the finite-element solution accuracies. These studies laid the foundation, or state of 
the art, for choosing optimum clement density in finite-element modeling of whole orbiter wing-type structures. 

Figure 6 shows the wing s lower skin temperature distributions predicted from the five thermal models. Except 
for model A, which had the same skin element density as that of the WING model (fig. 4), models B, C, D, and E 
gave very close values of the structural temperatures at the center region of the wing’s lower skin. 

Figure 7 shows how the number of radiation view factors F i; increases as the number of radiation elements 
R4 1 (or clement density) increases. For this particular type of wing cell structure, doubling the number of radiation 
elements R41 would cause the numberof radiation view factors F,y to nearly quadruple. 

Model E gives smooth structural temperature distribution (fig. 6), and its element density level can be considered 
an attractive element density to use in modeling the whole orbiter wing-type structures. However, the SPAR transient 
heat transfer analysis using model E required extremely long computer CPU time because of the computations of 
the time-dependent radiation exchange vector using the 93,868 radiation view factors F i} . Thus, the use of model 
E element density level in modeling the whole orbiter wing-type structures could be beyond the current computer 
capabilities (memory core space and computation speed). 

The structural temperature distribution given by model E has a fiat region near the wing skin central region in 
the chordwise direction. In this region, slightly coarser elements, as in the element density level of model D, could 
be used to reduce CPU time without sacrificing the solution accuracy. As shown in the Central Processing Unit 
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Time Optimization section of this report, the SPAR CPU time for model D is 21 to 27 percent of that for model 
E. The logical criteria for setting the clement density and size distribution for modeling large, complex aerospace 
structures could be somewhere between the element density levels of models D and E (for example, eight elements 
in chordwisc direction and six elements in spanwisc direction, with finer elements located near the heat sinks). 

CENTRAL PROCESSING UNIT TIME OPTIMIZATION 

In the finite-element transient heat transfer analysis using the SPAR program, the governing matrix equation 
used is of the form (ref. 1) 

(Kk+ K r + Kh)T + CT = Q + R+ H O) 


where 

Kk is the conduction matrix, 

K r the radiation matrix, 

Kh die convection matrix, 

T the absolute temperature, 

C the capacitance matrix, 

Q the source load vector, 

R the radiation load vector, 

H the convection load vector, and 

[ '] denotes time derivative. 

In the computations using equation (1), there are several key control parameters which could be adjusted to 
reduce the SPAR CPU time. However, at the same time the solution accuracy may be affected by adjusting those 
control parameters. Therefore, it is vital to find out the optimum values of those control parameters which will 
drastically reduce the SPAR CPU time and yet will hardly affect the solution accuracy. Those control param- 

ctcrs are: 

KTIME = control command to specify the time interval at which the system matrix K ( — Kk + K r + Kh) is to be 
factored (see eq. (1)). 

RTIME = control command to specify the time increments at which the radiation load vectors (R) are to be com- 
puted (see eq. (1)). 

RCONV = control command to set the convergence criterion for the radiation load vector computations. 

RNITER = control command to set the maximum number of iterations permitted during the computation of radia- 
tion load vectors (R). 

Figure 8 shows the SPAR transient thermal analysis (up to 3000 sec from reentry) computation CPU time (using 
the ELXSI 6400 computer) plotted as a function of the number of joint locations of thermal models for two sets of 
KTIME/RTIME values. By increasing both the values of KTIME and RTIME from 2 sec to 25 sec, the SPAR CPU 
time could be drastically reduced, and the solutions remain practically unchanged. The peak skin temperatures and 
actual SPAR CPU time associated with the different thermal models arc tabulated in table 1. Notice that the use o 
model E, which required a SPAR CPU time several times longer than that for model D, did not show any advantage 

in solution accuracy. 
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Tabic 1. SPAR computation CPU time and peak skin temperatures 
associated with different wing cell thermal models. 


Model 


KTIME 

= RTIME, sec 

SPAR CPU time, 
min 

Peak skin temperature, 
°F 

2 

25 

2 

25 

A 

19.38 

5.48 

101.11 

101.07 

B 

102.13 

23.18 

117.25 

117.27 

C 

177.87 

35.88 

117.17 

117.20 

D 

301.02 

56.45 

117.26 

117.30 

E 

1428.08 

204.28 

117.36 

117.42 


For conducting a more thorough study of the effects of KTIME, RTIME, and other SPAR control parameters 
previously mentioned, and for finding the optimum values of these parameters, the thermal model D (fig. 5) was 
used because it required much less SPAR CPU time than model E. 


Table 2. SPAR computation CPU times and peak skin temperatures based on thermal model D. 


KTIME, 

sec 

RTIME, 

sec 

RCONV 

RNITER 

SPAR CPU 
time, min 

Peak skin temperature, 
°F 

2 

2 

10 

0.0001 

12 

301.02 

117.26 

10 

231.08 

117.36 

25 

219.95 

117.44 

50 

218.15 

117.36 

100 

216.53 

116.55 

25° 

79.90 

117.39 

25“ 

56.45 

117.30 

10 

56.63 

117.30 

8 

56.90 

117.30 

5 6 

— 

C 

0.0005 

12 

53.98 

116.78 

0 . 001 6 

52.77 

116.24 

50 

0.0001 

50.10 

117.15 

100 

47.02 

116.87 

50 

50 

45.48 

117.07 


“Optimum KTIME and RTIME. 

& Default value. 

“Radiation exchange vectors did not converge. 


Table 2 summarizes the results of this investigation. Figure 9 shows the data of table 2 plotted in the KTIME/ 
RTIME space. Notice that by increasing KTIME from 2 sec to 25 sec (holding RTIME = 2 sec), the SPAR CPU time 
could be reduced from 301.02 min (point A) down to 219.95 min (point B). This is an approximately 27-percent 
reduction in SPAR CPU time. Further increase in the value of KTIME offered very little gain in saving SPAR CPU 
lime. Next, by increasing RTIME from 2 sec to 25 sec (keeping KTIME = 25 sec), the SPAR CPU time could be 
further reduced from 219.95 min (point B) down to 56.45 min (point Q. This gives another 54-percent reduction 
(based on point A) in the SPAR CPU time without affecting the solution accuracy. Further increase in the value of 
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RTIME beyond point C did not save much SPAR CPU time. At point C the SPAR CPU time is only 19 percent 
of that at point A, and the solution error is only 0.04 °F (compared to that at point A). Additional increase in both 
values of KTIME and RTIME to 50 sec (point D) had very little improvement in CPU time, and the solution error 
started to show up (0.19 °F error). Thus, point C could be chosen as the optimum point for selecting the values of 
KTIME and RTIME. The reductions of the values of RCONV and RNITER from 0.0001 and 12, respectively, had 
little effect on the improvement of the SPAR CPU time (table 2). 

MIDSPAN MODEL 

Element Density 

Based on the knowledge gained from the studies of element density requirements and the SPAR CPU time 
optimization described earlier, the element density criteria is established for the finite-element modelings of the 
orbiter wing-type structures. Before modeling the orbiter whole wing with high-density elements, it is wise to first 
model one segment of the orbiter wing with high-density elements and explore the CPU time requirements in the 
SPAR transient thermal analysis. The CPU time found for the wing segment model may then be extrapolated to 
estimate the CPU time requirement for the whole wing model having the same element density as the wing segment 
model. The said wing segment has four bays and is located at the midspan bounded by two adjacent wing ribs located 
at y o -226 and y o -254, respectively. 

The thermal model MIDSPAN for this wing segment is shown in figure 10. In the light of structural temperature 
distributions obtained from models D and E (fig. 6), the lower and upper skins of each bay of the wing segment 
were modeled with eight and six different-sized elements, respectively, in the chordwise and spanwise directions 
with finer elements used near the heat sinks (spar and rib caps). The MIDSPAN SPAR thermal model has 4064 joint 
locations, 1149 radiation elements R41, and 137,328 radiation view factors Fij. The corresponding NASTRAN 
(ref. 13) structural model has 722 grid points (18 percent of the thermal model). The sizes of the MIDSPAN thermal 
and structural models are compared in table 3. 


Table 3. Summary of thermal and structural models for 
an orbiter wing midspan segment. 


Item 

Model 

Thermal Structural 

Number of nodes 

4,046 

722 

Number of elements 

363 K21 
1,634 K41 
6 K61 
2,730 K81 
1,149 R41 
764 C41 

674 CQUAD2 
303 CROD 

Number of Fij 

137,328 

— 

Total CPU time, min 0 

464 

13 

Fij computation CPU time, min 

38 



°ELXSI 6400 computer 
KTIME = RTIME = 25 sec 
RCONV = 0.0001, RNITER = 12 
T 2 = 3000 sec 
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Central Processing Unit Time 


By using the previously established optimum values of KTIME = RTIME = 25 sec, transient heat transfer analysis 
was performed on the MIDSPAN thermal model. This thermal model required 464 min computation CPU time (table 
3) for the thermal analysis duration of 3000 sec. The corresponding structural model took 13 min computation CPU 
lime. It took only 38 min to compute the vast number (137,328) of radiation view factors. It must be emphasized 
that once the radiation view factors are computed, they are stored in the computer file, and no more computations 
are required in subsequent computer runs of the same thermal model. Therefore, in the transient thermal analysis 
most of the computation CPU time is used in computing time-dependent structural temperatures, not in computing 
radiation view factors, even though their number is enormous. This finding eliminated the fear of encountering a 
huge number of view factors in the thermal modelings of large, complex aerospace structures using finer elements. 

Structural Temperatures 

The calculations of the structural temperatures in the MIDSPAN thermal model were based on the STS-5 surface 
heating rates shown in figure 11. The distributions of skin tempc raturcs in time sequence are shown in figure 12. 
The most severe (highest values and gradients) temperature distributions occurred at 1700 sec from reentry, and 
gradually the severity tapered off as the lime increased. Judging from the temperature distribution surfaces for each 
bay, the use of an 8 by 6 system of element density gave acceptable smoothness of the skin temperature distributions. 
Therefore, an 8 by 6 system of element density could be considered as the element density criteria for modeling such 
types of aerospace structures. 

Thermal Stresses 

Based on the structural temperature distributions at 1700 sec, the thermal stresses ct x , a y , and r xy in the wing 
skins were calculated and are shown in figures 13 to 15, respectively. Again the distributions of all three stresses 
exhibited acceptable smoothness, emphasizing that 8 by 6 system element density criteria were satisfactory. Notice 
that a x , whose direction is transverse to the hat stringers of the wing skins, has a higher magnitude than that of <j y , 
and is therefore more critical. 


Thermal Deformations 

Figure 16 shows the undeformed and deformed shapes of the MIDSPAN structural model. The lower and upper 
skins of each bay bulged outwardly, with the maximum displacement of 0.08552 in. occurring near the midregion 
bay 1 ’s upper skin. This amount of deformation is way below the thickness (0.16 in.) of the SIP. A small region in 
bay l’s upper skin, near the spar cap separating bays 1 and 2, dented slightly inward. By observing figure 16, the 
MIDSPAN structural model gave an acceptable detailed deformation field. 

CONCLUSIONS 

The clement density level required in the finite-element thermal analysis of orbiler wing-lype structures was 
investigated, and the optimization of transient thermal analysis computation central processing unit (CPU) time was 
discussed. The newly established element density criteria for each wing cell were eight elements in the chordwise 
direction and six elements in the spanwisc direction. The resulting distributions of structural temperatures, thermal 
stresses, and thermal displacements were found to be satisfactory by using this level of clement density in the thermal 
modeling of an orbitcr wing segment. 
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The computation CPU time required for calculating the large number of radiation view factors comprised a 
small fraction of the total CPU time required for the computation of time-dependent structural temperatures. The 
study of the wing segment model, using the computer CPU time optimization technique, shows that the estimated 
compulation CPU time required for the thermal analysis of the whole orbitcr wing-type structures (with an 8 by 6 
element density system) may be within the acceptable range. 


Ames Research Center 
Dryden Flight Research Facility 
National Aeronautics and Space Administration 
Edwards, California, August 22 , 1988 
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Figure 1. Locations oflocal three-dimensional plug thermal models for space shuttle orbiter. 
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Figure 2 _ SPAR thermal model WS240 setup for orbiter wing midspan cross section. TPS elements removed to convert to structural model 



Thermal model 


Structural model 


Feature 

Number 

Feature 

Number 

JLOC 

605 

GRID point 

62 

K21 element 

250 

CQUAD2 element 

89 

K31 element 

6 

CROD element 

9 

K41 element 

516 



R21 element 

124 



C21 element 

96 



F.i 

1564 



* 






9468 


Figure 3. SPAR thermal model FS877 setup for orbiter midfuselage cross section. TPS and cargo bay door elements 
removed to convert to structural model (ref. 5). 


13 


■g CM CO O O) r- N 
t CO 0f> T" t— oo to 

3 CM t- 

Z 


c c c c c 

<D CD <D <D 4> 

EEEEE 

(D O <1) (t) (D 

Ooliio 

0 « lO r r Tt 

J W CM CO ^ 4 

-5 LU LU LU UJ LU 


CD 

•O O) CO ^ IT) 

P 00 O) CO 00 CO 

5 cm to 3F 

3 CMt- 


c c c c c 

CD 4) CD 4) 4> 

EEEEE 

CD CD 4) d) CD 

O 4> 4> 4) 4> 4) 

_J CM CO ** CO ** 

-5 * ^ * DC OC u. 



14 


Rgure 4 . SPAR thermal model WING setup for orbiter wing. TPS, wheel well door, and landing gear removed to convert to structural model 
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F U 

thermal 


K21 

K41 

K81 

R41 

C41 


model 








A 

112 

34 

28 

28 

15 

10 

78 

B 

436 

54 

168 

224 

89 

56 

2816 

C 

636 

82 

232 

336 

137 

88 

6894 

D 

972 

98 

360 

560 

201 

120 

13500 

E 

2076 

146 

848 

1344 

513 

320 

93869 

NASTRAN 

Grid 


CQUAD2 

CROD 




structural 








model 








A 

24 


18 

54 




B 

82 
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Figure 5. SPAR thermal models A, B, C, D, E setups for orbiter wing midspan bay 3. The K81 elements for TPS 
and SIP not shown; TPS and SIP removed to convert to structural models. 
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Figure 6. Distributions of orbiter wing lower skin temperatures predicted from thermal models A, B, C, D, E; time 
1,700 sec, STS-5. 
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Figure 7. Number of radiation view factors as a function of number of radiation elements, orbiter wing midspan 
3. 



KTIME = RTIME = 2 sec 
KTIME = RTIME = 25 sec 
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Figure 9. Three-dimensional plots of SPAR CPU time in KTIME— RTIME space. 
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Figure 10. SPAR thermal model MIDSPAN setup for orbiter wing midspan segment; TPS removed to convert to structural model. 
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Figure 12. Orbiter wing skin temperature distributions 
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Figure 12. Continued. 
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Figure 12. Continued. 
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Time = 2000 sec 

Figure 12. Continued. 
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Figure 12. Continued. 
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Figure 12. Concluded. 







Figure 13. Distributions of chordwisc stress a x in orbiter wing skins; time = 1700 sec, STS-5. 
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Figure 15. Distributions of shear stress r xy in orbiter wing skins; time = 1700 sec, STS-5. 
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Figure 16. Deformed shape of orbiter wing midspan segment; time = 1700 sec, STS-5. 
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